*FIG 5: EVENT STUDIES ACROSS ALL CANTONS *

// NOTE: data needed for this do-file cannot be shared due to privacy reasons

/* following packages need to be installed:
ssc install winsor2
ssc install xtvent
ssc install coefplot
*/

version 16.1

clear all
cap log close 
cap clear matrix
set more off
cap set scheme mygraphs


* * *   * * *   * * *   * * *   * * *   * * *
*   LOAD & PREPARE THE DATA   
* * *   * * *   * * *   * * *   * * *   * * *
run "$mypathRR/Resources/prep_municipality_level_data.do"

drop if (year == 1996 | year == 1998 | year == 1999 | year == 2000)

egen periods = group(year)

xtset gemeinde periods

gen policyvar = (cant == 6 & year>2005)

global outcomes "share300K share300Kc share300Kc2 reink reink_tr reink_ctr reink_w reink_tr_w reink_ctr_w share200K  " 
global outcomespaper "share300K  reink_w  share200K " 


// weighted regression, cluster canton-year and gemeinde, 2 models in 1 plot
foreach outcome in $outcomespaper {

	* model 1: pre-trend correction
	xtevent `outcome' [aweight = weight_`outcome'], policyvar(policyvar) 	///
	window(-7 9) trend(-7) reghdfe  										///
	vce(cluster cant#periods gemeinde)  
	est sto `outcome'_pretrend_II
	
	xteventplot, name(comp_`outcome'II, replace) nosupt minus1label 		///
	ytitle("`: variable label `outcome''") 									///
	ylabel(`myla', labsize(small) grid) 									///
	xlab( , labsize(small) ang(45)) 										///
	xline(-0.2, lcolor(red)) xline(1.8, lcolor(green)) 

	* model 2: no pre-trend correction
	xtevent `outcome' [aweight = weight_`outcome'], policyvar(policyvar) 	///
	window(-7 9) reghdfe  													///
	vce(cluster cant#periods gemeinde)  
	est sto `outcome'
	
	xteventplot, name(comp_`outcome', replace) nosupt minus1label 			///
	ytitle("`: variable label `outcome''") 									///
	ylabel(`myla', labsize(small) grid) 									///
	xlab(, labsize(small) ang(45)) 											///
	xline(-0.2, lcolor(red)) xline(1.8, lcolor(green)) 
}

* * * * PLOT GRAPHS * * * * 

global coeflistxtevent " _k_eq_m8 _k_eq_m7 _k_eq_m6 _k_eq_m5 _k_eq_m4 _k_eq_m3 _k_eq_m2  _k_eq_p0 _k_eq_p1 _k_eq_p2 _k_eq_p3 _k_eq_p4 _k_eq_p5 _k_eq_p6 _k_eq_p7 _k_eq_p8 _k_eq_p9 _k_eq_p10"
global coeflistxteventtrend "   _k_eq_p0 _k_eq_p1 _k_eq_p2 _k_eq_p3 _k_eq_p4 _k_eq_p5 _k_eq_p6 _k_eq_p7 _k_eq_p8 _k_eq_p9 _k_eq_p10"

/* 	estimates stored in e(b) OR e(delta)
	covariance matrix is e(V) OR e(Vdelta)	
*/
foreach OUT in $outcomespaper {
		estimates restore `OUT'
		
		* manually calculate confidence intervals from stored results
		loc df = e(df)
		di `df'
		tempname b V
		mat `b' = e(delta)
		mat `V' = e(Vdelta)
		tempvar se	
		
		* standard erors
	    mata st_matrix("`se'", sqrt( diagonal(st_matrix("`V'"))) ')	

		mat list `b'
		mat list `se'
		
		* confidence intervals
		forval n = 1/18 {
		loc ll`n' = `b'[1,`n'] - invttail(`df',0.5*0.05) * `se'[1,`n'] 
		loc ul`n' = `b'[1,`n'] + invttail(`df',0.5*0.05) * `se'[1,`n'] 
		}

		* create a new matrix to put the results
		matrix A`OUT' = J(3,19,.)
        matrix rownames A`OUT' = est ll95 ul95
        matrix colnames A`OUT' = _k_eq_m8 _k_eq_m7 _k_eq_m6 _k_eq_m5 _k_eq_m4 _k_eq_m3 _k_eq_m2 _k_eq_m1 _k_eq_p0 _k_eq_p1 _k_eq_p2 _k_eq_p3 _k_eq_p4 _k_eq_p5 _k_eq_p6 _k_eq_p7 _k_eq_p8 _k_eq_p9 _k_eq_p10
		mat list A`OUT'
		
		* put the results from xtevent into the new matrix
		forval n = 1/7 {
		loc i = `n'
		matrix A`OUT'[1,`n'] = `b'[1,`i'] \ `ll`i'' \ `ul`i''
		}
	 	matrix A`OUT'[1,8] = 0 \ 0 \ 0
		forval n = 9/19 {
		loc i = `n'-1
		matrix A`OUT'[1,`n'] = `b'[1,`i'] \ `ll`i'' \ `ul`i''
		}
		
		* compare results stored in new matrix with output from xtevent
		mat list A`OUT'
		mat list `se'
		estimates replay `OUT'
	 }
	 
foreach OUT in $outcomespaper {
		estimates restore `OUT'_pretrend_II
		
		* manually calculate confidence intervals from stored results
		loc df = e(df)
		di `df'
		tempname b V
		mat `b' = e(delta)
		mat `V' = e(Vdelta)
		tempvar se	
		
		* standard erors
	    mata st_matrix("`se'", sqrt( diagonal(st_matrix("`V'"))) ')	

		mat list `b'
		mat list `se'
		
		* confidence intervals
		forval n = 1/12 {
		loc ll`n' = `b'[1,`n'] - invttail(`df',0.5*0.05) * `se'[1,`n'] 
		loc ul`n' = `b'[1,`n'] + invttail(`df',0.5*0.05) * `se'[1,`n'] 
		}

		* create a new matrix to put the results
		matrix B`OUT' = J(3,19,.)
        matrix rownames B`OUT' = est ll95 ul95
        matrix colnames B`OUT' = _k_eq_m8 _k_eq_m7 _k_eq_m6 _k_eq_m5 _k_eq_m4 _k_eq_m3 _k_eq_m2 _k_eq_m1 _k_eq_p0 _k_eq_p1 _k_eq_p2 _k_eq_p3 _k_eq_p4 _k_eq_p5 _k_eq_p6 _k_eq_p7 _k_eq_p8 _k_eq_p9 _k_eq_p10
		mat list B`OUT'
		
		* put the results from xtevent into the new matrix
		forval n = 1/1 {
		loc i = `n'
		matrix B`OUT'[1,`n'] = `b'[1,`i'] \ `ll`i'' \ `ul`i''
		}
		forval n = 2/8 {
		matrix B`OUT'[1,`n'] = 0 \ 0 \ 0
		}
		forval n = 9/19 {
		loc i = `n'-7
		matrix B`OUT'[1,`n'] = `b'[1,`i'] \ `ll`i'' \ `ul`i''
		}
		
		* compare results stored in new matrix with output from xtevent
		mat list B`OUT'
		mat list `se'
		estimates replay `OUT'_pretrend_II
	 }
	  
	 
// graphs with coefficients in percentage changes
preserve 
foreach outcome in $outcomespaper {
	
	* gen avg. in 2005 for relative changes
	sum `outcome' [aweight = weight_`outcome'] if cant == 6 & year == 2005
	global `outcome'2005 = `r(mean)'
	
	* outcome-specific labels 
	if "`outcome'" == "share300K" | "`outcome'" == "share300Kc"| "`outcome'" == "share300Kc2"| "`outcome'" == "share200K"  {
		label var `outcome' "Rel. change in share of top earners (in %)"
		local mylabel "0(25)80"
	} 
	if "`outcome'" == "share200K"   {
		label var `outcome' "Rel. change in share of taxpayers 200K<y<300K (in %)"
	} 
	 if "`outcome'" == "reink" | "`outcome'" == "reink_tr" | "`outcome'" == "reink_ctr" | "`outcome'" == "reink_w" | "`outcome'" == "reink_tr_w" | "`outcome'" == "reink_ctr_w" {
		label var `outcome' "Change in net income per taxpayer (in %)"
		local mylabel "-15(5)70"
	}
	else {
		local mylabel
	}

	* plot coefficients
		coefplot    (matrix(A`outcome'), ci((2 3)) label(" ") 					///
					omitted color(*0.2) ciopts(color(*0.3)) connect(direct)		///
					transform(* = @/${`outcome'2005}*100 ))						///
					(matrix(B`outcome'), ci((2 3)) label(" ") omitted connect(direct)			///
					 transform(* = @/${`outcome'2005}*100 )) 									///
					, coeflabels( _k_eq_m8 = "-1993"  _k_eq_m7 = "1995" _k_eq_m6 = "1997" 		///
					_k_eq_m5  = "2001" _k_eq_m4 = "2002" _k_eq_m3 = "2003" _k_eq_m2 = "2004"	///  
					_k_eq_m1 = "2005" 															///
					_k_eq_p0 = "2006" _k_eq_p1 = "2007" _k_eq_p2 = "2008" _k_eq_p3 = "2009" 	///
					_k_eq_p4 = "2010" _k_eq_p5 = "2011" _k_eq_p6 = "2012" _k_eq_p7 = "2013" 	///
					_k_eq_p8 = "2014" _k_eq_p9 = "2015"  _k_eq_p10 = "2016") 					///
					vertical levels(95) ciopts(recast(rcap))  												///
					yline(0, lcolor(black))  xline(8.7, lcolor(red)) xline(10.7, lcolor(green)) 				///
					legend(order(2 "w/o controlling for pre-trend" 4 "controlling for pre-trend") 			///
					  row(2) ring(0) position(11) bmargin(small) size(small) ) 								///
					xlabel(, labsize(small) ang(45)) 						///
					ylabel(`mylabel', labsize(small) grid) 					///
					ytitle("`: variable label `outcome''") 					///
					name(gr_`outcome'_pct, replace) 
					
		if "`outcome'" == "reink_w" {
		graph export "Results/Fig_5b)-`outcome'_0_wght_cluster-pct_change-xtevent.pdf", as(pdf) replace
		}
											
}
restore



// graphs with coefficients in levels
preserve 
foreach outcome in $outcomespaper {
	
	* gen avg. in 2005 
	sum `outcome' [aweight = weight_`outcome'] if cant == 6 & year == 2005
	global `outcome'2005 = `r(mean)'
	local mean : di %3.2f `=${`outcome'2005}'
	
	* outcome-specific labels 
	if "`outcome'" == "share300K" | "`outcome'" == "share300Kc"| "`outcome'" == "share300Kc2"| "`outcome'" == "share200K"  {
		label var `outcome' "Effect on share of top earners (in pp)"
		local mylabel "0(0.1)1"
		local textplace "0.75 3"
	} 
	if "`outcome'" == "share200K"   {
		label var `outcome' "Effect on share of taxpayers with 200K < y < 300K (in pp)" 
		local mylabel "0(0.1)1"
		local textplace "0.75 3"

	} 
	 if "`outcome'" == "reink" | "`outcome'" == "reink_tr" | "`outcome'" == "reink_ctr" | "`outcome'" == "reink_w" | "`outcome'" == "reink_tr_w" | "`outcome'" == "reink_ctr_w" {
		label var `outcome' "Effect on net income per taxpayer (in CHF)"
		local mylabel "0(5)40"
		local textplace "32 3"

	}
	
	* plot coefficients
		coefplot    (matrix(A`outcome'), ci((2 3)) label(" ") 					///
					omitted color(*0.2) ciopts(color(*0.3)) connect(direct)	)	///
					(matrix(B`outcome'), ci((2 3)) label(" ") omitted connect(direct)	)	 	///
					, coeflabels( _k_eq_m8 = "-1993"  _k_eq_m7 = "1995" _k_eq_m6 = "1997" 		///
					_k_eq_m5  = "2001" _k_eq_m4 = "2002" _k_eq_m3 = "2003" _k_eq_m2 = "2004"	///  
					_k_eq_m1 = "2005" 															///
					_k_eq_p0 = "2006" _k_eq_p1 = "2007" _k_eq_p2 = "2008" _k_eq_p3 = "2009" 	///
					_k_eq_p4 = "2010" _k_eq_p5 = "2011" _k_eq_p6 = "2012" _k_eq_p7 = "2013" 	///
					_k_eq_p8 = "2014" _k_eq_p9 = "2015"  _k_eq_p10 = "2016") 					///
					vertical levels(95) ciopts(recast(rcap))  												///
					yline(0, lcolor(black))  xline(8.7, lcolor(red)) xline(10.7, lcolor(green)) 			///
					legend(order(2 "w/o controlling for pre-trend" 4 "controlling for pre-trend") 			///
					row(2) ring(0) position(11) bmargin(small) size(small) ) 								///
					text(`textplace' "avg. 2005 (OW) = `mean'", size(small))								///
					xlabel(, labsize(small) ang(45)) 						///
					ylabel(`mylabel', labsize(small) grid) 					///
					ytitle("`: variable label `outcome''") 					///
					name(gr_`outcome'_pct, replace)
		
		if "`outcome'" == "share300K" {
		graph export "Results/Fig_5a)-`outcome'_0_wght_cluster-xtevent.pdf", as(pdf) replace
		}
}

restore


						* * * * *  E N D  * * * * * * 
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
